##########
# R Info #
##########
version
#platform       x86_64-w64-mingw32          
#arch           x86_64                      
#os             mingw32                     
#system         x86_64, mingw32             
#status                                     
#major          3                           
#minor          6.0                         
#year           2019                        
#month          04                          
#day            26                          
#svn rev        76424                       
#language       R                           
#version.string R version 3.6.0 (2019-04-26)
#nickname       Planting of a Tree 

library(ggplot2)
library(ggpubr)
library(maps)
library(wesanderson)
library(stargazer)
library(dplyr)

sessionInfo()
#R version 3.6.0 (2019-04-26)
#Platform: x86_64-w64-mingw32/x64 (64-bit)
#Running under: Windows 10 x64 (build 19041)
#Matrix products: default
#locale:
#[1] LC_COLLATE=English_United States.1252  LC_CTYPE=English_United States.1252   
#[3] LC_MONETARY=English_United States.1252 LC_NUMERIC=C                          
#[5] LC_TIME=English_United States.1252    
#attached base packages:
#[1] stats     graphics  grDevices utils     datasets  methods   base     
#other attached packages:
#[1] dplyr_1.0.4       stargazer_5.2.2   wesanderson_0.3.6 maps_3.3.0        ggpubr_0.4.0     
#[6] ggplot2_3.3.2    
#loaded via a namespace (and not attached):
#[1] zip_2.0.4         Rcpp_1.0.3        cellranger_1.1.0  pillar_1.4.3      compiler_3.6.0   
#[6] forcats_0.5.0     tools_3.6.0       lifecycle_0.2.0   tibble_3.0.4      gtable_0.3.0     
#[11] pkgconfig_2.0.3   rlang_0.4.10      openxlsx_4.1.4    DBI_1.1.0         rstudioapi_0.11  
#[16] curl_4.3          haven_2.2.0       rio_0.5.16        rJava_0.9-11      withr_2.4.1      
#[21] dplyr_1.0.4       hms_0.5.3         generics_0.0.2    xlsxjars_0.6.1    vctrs_0.3.6      
#[26] grid_3.6.0        tidyselect_1.1.0  data.table_1.12.8 glue_1.4.1        R6_2.4.1         
#[31] rstatix_0.6.0     readxl_1.3.1      foreign_0.8-71    carData_3.0-3     car_3.0-7        
#[36] purrr_0.3.3       tidyr_1.1.2       magrittr_1.5      backports_1.1.5   scales_1.1.0     
#[41] ellipsis_0.3.0    abind_1.4-5       assertthat_0.2.1  colorspace_1.4-1  ggsignif_0.6.0   
#[46] xlsx_0.6.1        stringi_1.4.5     munsell_0.5.0     broom_0.7.5       crayon_1.3.4




#############
# Load Data #
#############
load("PostDataJQD.RData")

colnames(df)
# country:          country name
# year:             year of posting
# yearmonth:        year + month of posting
# party_country:    party name + country name
# party:            party name
# cmp_familyid:     party family id in the CMP
# cmp_family:       party family name in the CMP
# ches_ideology:    general left-right ideology in the CHES
# gps_populism:     populism score in the GPS
# gps_econlr:       economic left-right ideology in the GPS      
# gps_soclr:        social left-right ideology in the GPS
# Page.Name:        name of facebook page
# Likes.at.Posting: number of page likes at posting
# Likes:            number of Likes
# Comments:         number of comments
# Shares:           number of shares
# Love:             number of Love reactions
# Wow:              number of Wow reactions
# Haha:             number of Haha reactions 
# Sad:              number of Sad reactions
# Angry:            number of Angry reactions
# totalreactions:   sum of Likes, Love, Wow, Haha, Sad, and Angry
# loveprop:         proportion of Love reactions (Love / totalreactions)
# angryprop:        proportion of Angry reactions (Angry / totalreactions)

summary(df)

# number of countries
length(unique(df$country)) # 79

# number of parties
length(unique(df$party_country)) # 690

# compute emotional polarization
df$lapolarization <- (df$angryprop + 1) / (df$loveprop + 1)

# set colors
rgy <- wes_palette("Darjeeling1", 3)




############################################
# Figure 2: Love and Angry Country Ranking #
############################################
# compute average love proportion, angry proportion, and emotional polarization by country
countries <- data.frame(country=sort(unique(df$country)))
countries$aveloveprop <- tapply(df$loveprop, df$country, mean)
countries$aveangryprop <- tapply(df$angryprop, df$country, mean)
countries$lapolarization <- tapply(df$lapolarization, df$country, mean)

# function to get confidence intervals
getCI <- function(vector){
  m <- mean(vector)
  s <- sd(vector)
  error <- qt(0.995,df=length(vector)-1) * s / sqrt(length(vector))
  left <- m - error
  right <- m + error
  return(c(left, right))
}

#pdf("AverageLoveAngry.pdf", width=12, height=15)
par(mfrow=c(1,2), mar=c(3.1,8,3.1,2.1))

# order countries by average love proportion
countries <- countries[order(countries$aveloveprop),]

# figure 2 left panel
plot(countries$aveloveprop, 1:79, type="n",
     ylab="", yaxt="n", ylim=c(3,77),
     xlab="", xlim=c(0,0.13),
     main="Average Love Proportion")
abline(h=1:79, col="gray95")
abline(v=mean(df$loveprop), lty=2)
points(countries$aveloveprop, 1:79, pch=20, cex=2, col=rgy[2])
axis(2, 1:79, labels=countries$country, las=1, cex.axis=0.8)
text(mean(df$loveprop),1, labels="Global Mean", cex=0.7)
for(i in 1:79){
  ci <- getCI(df$loveprop[df$country == countries$country[i]])
  segments(ci[1], i, ci[2], i, lwd=1.3, col=rgy[2])
}


# order countries by average love proportion
countries <- countries[order(countries$aveangryprop),]

# figure 2 right panel
plot(countries$aveangryprop, 1:79, type="n",
     ylab="", yaxt="n", ylim=c(3,77),
     xlab="", xlim=c(0,0.13),
     main="Average Angry Proportion")
abline(h=1:79, col="gray95")
abline(v=mean(df$angryprop), lty=2)
points(countries$aveangryprop, 1:79, pch=20, cex=2, col=rgy[1])
axis(2, 1:79, labels=countries$country, las=1, cex.axis=0.8)
text(mean(df$angryprop),1, labels="Global Mean", cex=0.7)
for(i in 1:79){
  ci <- getCI(df$angryprop[df$country == countries$country[i]])
  segments(ci[1], i, ci[2], i, lwd=1.3, col=rgy[1])
}

#dev.off()




####################################################
# Figure 3: Emotional Polarization Country Ranking #
####################################################
# order countries by average emotional polarization
countries <- countries[order(countries$lapolarization),]

#pdf("PolarizationRank.pdf", width=7, height=15)
par(mfrow=c(1,1), mar=c(3.1,9,3.1,2.1))

# figure 3
plot(countries$lapolarization, 1:79, type="n",
     ylab="", yaxt="n", ylim=c(3,77),
     xlab="", xlim=c(0.9,1.1),
     main="Average Emotional Polarization")
abline(h=1:79, col="gray95")
abline(v=1, lty=2)
points(countries$lapolarization, 1:79, pch=20, cex=2, col=rgy[3])
axis(2, 1:79, labels=countries$country, las=1, cex.axis=0.8)
for(i in 1:79){
  ci <- getCI(df$lapolarization[df$country == countries$country[i]])
  segments(ci[1], i, ci[2], i, lwd=1.3, col=rgy[3])
}

#dev.off()




##################
# Figure 4: Maps #
##################
world_map <- map_data("world")

# check countries do not match
countries$country[which(!(countries$country %in% unique(world_map$region)))]
# United Kingdom Bosnia Herzegovina North Macedonia United States 

countries$region <- as.character(countries$country)

# fix country names
countries$region[countries$country == "Bosnia Herzegovina"] <- "Bosnia and Herzegovina"
countries$region[countries$country == "North Macedonia"] <- "Macedonia"
countries$region[countries$country == "United Kingdom"] <- "UK"
countries$region[countries$country == "United States"] <- "USA"

which(!(countries$region %in% world_map$region)) # integer(0)

# subset data
loveangry.map <- right_join(countries[,c("region", "aveloveprop", "aveangryprop", 
                                         "lapolarization")], 
                            world_map, by = "region")

# exclude antarctica
loveangry.map <- loveangry.map[loveangry.map$region != "Antarctica",]

# change column names
colnames(loveangry.map)[2:4] <- c("Ave Love Prop",
                                  "Ave Angry Prop",
                                  "Ave EP")

# group average emotional polarization by bin
loveangry.map$`Ave EP` <- cut(loveangry.map$`Ave EP`,
                              c(1.100,1.050,1.025,1.000,
                                0.975,0.950,0.900), dig.lab=4)
loveangry.map$`Ave EP` <- factor(loveangry.map$`Ave EP`,
                                 c("(1.05,1.1]", "(1.025,1.05]", 
                                   "(1,1.025]", "(0.975,1]", 
                                   "(0.95,0.975]", "(0.9,0.95]"))

# set ggplot themes
plain <- theme(
  axis.text = element_blank(),
  axis.line = element_blank(),
  axis.ticks = element_blank(),
  panel.border = element_blank(),
  panel.grid = element_blank(),
  axis.title = element_blank(),
  panel.background = element_rect(fill = "white"),
  plot.title = element_text(hjust = 0.5)
)

plain2 <- theme(
  axis.text = element_blank(),
  axis.line = element_blank(),
  axis.ticks = element_blank(),
  panel.border = element_rect(fill = NA),
  panel.grid = element_blank(),
  axis.title = element_blank(),
  panel.background = element_rect(fill = "white"),
  plot.title = element_text(hjust = 0.5),
  legend.position = "none"
)

# figure 4 top panel
worldlove <- ggplot(data = loveangry.map, mapping = aes(x = long, y = lat, group = group)) + 
  geom_polygon(aes(fill = `Ave Love Prop`), color = "white", size = 0.1) +
  scale_fill_distiller(palette ="Greens", direction = 1, limits = c(0, 0.13)) + 
  geom_rect(xmin = -13.8, xmax = 32.5, ymin = 33, ymax = 74, 
            fill = NA, colour = "black", size = 0.8) +
  plain

europelove <- ggplot(data = loveangry.map, mapping = aes(x = long, y = lat, group = group)) + 
  geom_polygon(aes(fill = `Ave Love Prop`), color = "white", size = 0.1) +
  scale_fill_distiller(palette ="Greens", direction = 1, limits = c(0, 0.13)) +
  coord_sf(xlim = c(-13, 32.5), ylim = c(33, 74), expand = FALSE) +
  plain2

lovemp <- ggplot() +
  coord_equal(xlim = c(0, 3.3), ylim = c(0, 1), expand = FALSE) +
  annotation_custom(ggplotGrob(worldlove), xmin = 0, xmax = 2.3, ymin = 0, 
                    ymax = 1) +
  annotation_custom(ggplotGrob(europelove), xmin = 2.3, xmax = 3, ymin = 0, 
                    ymax = 1) +
  ggtitle(label = "A. Average Love Proportion") +
  theme_void() +
  theme(plot.title = element_text(hjust = 0.5, face = "bold"))


# figure 4 middle panel
worldangry <- ggplot(data = loveangry.map, mapping = aes(x = long, y = lat, group = group)) + 
  geom_polygon(aes(fill = `Ave Angry Prop`), color = "white", size = 0.1) +
  scale_fill_distiller(palette ="Reds", direction = 1, limits = c(0, 0.13)) + 
  geom_rect(xmin = -13.8, xmax = 32.5, ymin = 33, ymax = 74, 
            fill = NA, colour = "black", size = 0.8) +
  plain

europeangry <- ggplot(data = loveangry.map, mapping = aes(x = long, y = lat, group = group)) + 
  geom_polygon(aes(fill = `Ave Angry Prop`), color = "white", size = 0.1) +
  scale_fill_distiller(palette ="Reds", direction = 1, limits = c(0, 0.13)) +
  coord_sf(xlim = c(-13, 32.5), ylim = c(33, 74), expand = FALSE) +
  plain2

angrymp <- ggplot() +
  coord_equal(xlim = c(0, 3.3), ylim = c(0, 1), expand = FALSE) +
  annotation_custom(ggplotGrob(worldangry), xmin = 0, xmax = 2.3, ymin = 0, 
                    ymax = 1) +
  annotation_custom(ggplotGrob(europeangry), xmin = 2.3, xmax = 3, ymin = 0, 
                    ymax = 1) +
  ggtitle(label = "B. Average Angry Proportion") +
  theme_void() +
  theme(plot.title = element_text(hjust = 0.5, face = "bold"))


# figure 4 bottom panel
worlddiff2 <- ggplot(data = loveangry.map, mapping = aes(x = long, y = lat, group = group)) +
  geom_polygon(aes(fill = `Ave EP`), color = "white", size = 0.1) +
  scale_fill_brewer(palette ="RdYlGn", direction = 1, na.value = "grey50") + 
  geom_rect(xmin = -13.8, xmax = 32.5, ymin = 33, ymax = 74, 
            fill = NA, colour = "black", size = 0.8) +
  plain

europediff2 <- ggplot(data = loveangry.map, mapping = aes(x = long, y = lat, group = group)) + 
  geom_polygon(aes(fill = `Ave EP`), color = "white", size = 0.1) +
  scale_fill_brewer(palette ="RdYlGn", direction = 1, na.value = "grey50") + 
  coord_sf(xlim = c(-13, 32.5), ylim = c(33, 74), expand = FALSE) +
  plain2

diffmp2 <- ggplot() +
  coord_equal(xlim = c(0, 3.3), ylim = c(0, 1), expand = FALSE) +
  annotation_custom(ggplotGrob(worlddiff2), xmin = 0, xmax = 2.3, ymin = 0, 
                    ymax = 1) +
  annotation_custom(ggplotGrob(europediff2), xmin = 2.3, xmax = 3, ymin = 0, 
                    ymax = 1) +
  ggtitle(label = "C. Average Emotional Polarization") +
  theme_void() +
  theme(plot.title = element_text(hjust = 0.5, face = "bold"))


# figure 4 all
ggarrange(lovemp, angrymp, diffmp2,
          ncol = 1, nrow = 3)
#ggsave("LoveAngryMap.pdf", width=14, height=14)
